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One essential ingredient of evolutionary theory is the concept of htness as a measure for a species’ 
success in its living conditions. Here, we quantify the effect of environmental fluctuations onto fitness 
by analytical calculations on a general evolutionary model and by studying corresponding individual- 
based microscopic models. We demonstrate that not only larger growth rates and viabilities, but 
also reduced sensitivity to environmental variability substantially increases the fitness. Even for 
neutral evolution, variability in the growth rates plays the crucial role of strongly reducing the 
expected fixation times. Thereby, environmental fluctuations constitute a mechanism to account 
for the effective population sizes inferred from genetic data that often are much smaller than the 
census population size. 
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Spencer’s famous expression “survival of the fittest” [T] 
provides an appealing short summary of Darwin’s con¬ 
cept of evolution ElE]. However, it leaves aside a very 
difficult yet important aspect namely identifying the fac¬ 
tors determining the fitness of a species [H [S] : fittest 
individuals are by definition prevailing but the reasons 
facilitating their survival are not obvious. Besides the 
difficulties arising due to the genotype phenotype map¬ 
ping causing complex fitness landscapes [6], also ecologi¬ 
cal factors like population structure and composition ad¬ 
ditionally complicate the issue. Therefore traditional fit¬ 
ness concepts solely based on growth rates and viability 
were extended by frequency-dependent or inclusive 
fitness approaches [8]. Another important factor for the 
success of a certain trait, is a non-constant environment 
influencing birth/death rates [51418). How variable en¬ 
vironmental conditions affect evolutionary strategies like 
phenotypic heterogeneity or bet-hedging has been exten¬ 
sively studied, see e.g. [IS1IIM23], yet the consequences 
of fluctuating reproduction rates and their interplay with 
demographic fluctuations were not fully elucidated. 

Here, we quantitatively investigate the impact of vari¬ 
able environments on the fitness. In contrast to other 
models dealing with variable environmental conditions, 
we do not study which strategy is optimally suited to 
cope with such changing environments but focus on the 
consequences of fluctuating reproduction rates. In par¬ 
ticular, we show that an individual’s sensitivity to en¬ 
vironmental changes contributes substantially to its fit¬ 
ness: A reduced sensitivity increases the htness and may 
compensate for large disadvantages in the average repro¬ 
duction rate. We also hnd that huctuating environments 
inhuence neutral evolution where they can cause much 
quicker hxation times than expected. These effects are 
relevant as constant environmental conditions are the ex¬ 
ception rather than the norm; for instance, the availabil¬ 
ity of different nutrients, the presence of detrimental sub¬ 
stances and other external factors like temperature, all 
strongly inhuence reproduction/survival and occur on a 


broad range of time scales [51] . 

To understand the impact of huctuating environmen¬ 
tal conditions, we hrst consider rapidly changing envi¬ 
ronments in an evolutionary process based on birth and 
death events similar to [12]. The dynamics is described 
by the following stochastic differential equations: 


Ns = \vs-1-^\Ns+Ns(js^s + \INs{vs+i-^)us- ( 1 ) 

The inhuence of environmental variability is modeled as 
white noise acting on the growth rate, ug, of a trait of 
type S-. vs + where = 5{t - t') and 

as is the Standard Deviation (STD) of the noise. Death 
rates are assumed to be constant and identical for all 
traits [25|. Population growth is bounded and there¬ 
fore death rates increase with the total population size, 
N = '^s where Ns is the number of 5'-type indi¬ 
viduals [26j . This may account for density-dependent 
ecological factors such as limited resources or metabolic 
waste products accumulating at high population sizes. 
For specihcity, we choose 'yN/K as functional form where 
K is the carrying capacity scaling the maximal number of 
individuals and 7 sets the rate of death events. Beside en¬ 
vironmental noise, demographic huctuations arising from 
the stochastic nature of the birth-death dynamics yield 
the term y/Ns{vs+ uN/K)tis, where fis is (^-correlated 
noise, {pisit)Psit')) =S{t — t'), with a variance given by 
the sum of reaction rates m- Both multiplicative noise 
terms in Eq. Q are interpreted in the Ito sense [2H|. Note 
that environmental noise is linearly multiplicative in Ns, 
which is crucial for our results. 

Let us now consider the Fokker-Planck equation (FPE) 
associated to Eq. 0 , which we will use to derive fixation 
probabilities and times. We shall carry out further anal¬ 
ysis for two different traits S G {I52}; generalizations 
are straightforward. The transformation of Eq. ([^ to 
a FPE, depends on the correlation level of the environ¬ 
mental noise acting on distinct traits [29]. While demo¬ 
graphic noise for different traits is always uncorrelated. 
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the same environmental noise can affect multiple traits, 
e.g. if both traits feed from the same nutrients whose 
abundance fluctuates. We keep the following analysis 
quite general by introducing the correlation coefficient e : 
for e = 0 environmental noise is uncorrelated, (C 1 C 2 ) = 0 , 
while e = ±l for ^i=±^ 2 - The resulting FPE is 


dPiNi,N2,t) 

dt 


= edl 2(^i(T2NiN2P 


N 


+ 







( 2 ) 
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selection strength, s 


where di = . To uncover the influence of environ¬ 

mental noise on the evolutionary dynamics, the relative 
abundances seem the natural observables. Therefore, we 
change variables to the fraction x = and the to¬ 

tal number of individuals N = Ni + N 2 . The FPE for 
X and N can be simplified exploiting the fact that se¬ 
lection, s = vi — V 2 ^ is much slower than population 
growth vix 122(1 — x). Therefore, we integrate over 
the total population size N, considering the FPE for 
the marginal distribution P{x) = P(x, N) dN, and 
employing iV 3> 1, see SM [30]. The resulting one¬ 
dimensional FPE reads: 


dP{x, t) 
Wt 



dx { [—s —cr2(l~2;)-|-crJa::-|-eCTicr2 (1 — 2x)] Q} 


-2ecricr2+c’'2 




Q\=CP{x,t), (3) 


where Q = a;(l — x)P{x,t) and the last equality defines 
the Fokker-Planck operator C needed in the sequel. In 
Ref. [To] a similar FPE was derived for the special case 
(Ti = 02 - For 01—02 = 0, the drift term reduces to 
the well-know expression —sdxQ favoring the trait with 
a higher growth rate [5|. Note that the variability in 
the growth rates affects both the diffusion term and the 
drift, which is due to the multiplicative nature of the 
environmental noise. 

For simplicity, we discuss the case of different envi¬ 
ronmental sensitivity, defined by cti = A and 02 = 0, 
i.e. only the reproduction rate of the first trait depends 
on the environment. The drift is then proportional to 
a{x) = (s — A^a:)a:(l — x) independent of e. If s < 0, i.e. 
the second trait with a smaller variability in its birth rate 
is also faster in reproducing, the evolutionary dynamics 
does not change qualitatively compared to A = 0. Con¬ 
versely, if s > 0 , the situation changes dramatically: the 
growth rate favors trait 1 while the variability term fa¬ 
vors trait 2 . This leads to a stable fixed point x* = -^ 
for s < A^ (for s > A^ variability is not sufficient to pre¬ 
vent extinction of trait 2 ). Such a dynamics can be inter¬ 
preted as frequency-dependent fitness function. However, 
the frequency-dependence arises here from environmental 
noise and not from a pay-off matrix [32] as in standard 
evolutionary game theory. 


FIG. 1: Fixation probability, Pax, depending on selection 
strength, s, and variability cri = A according to Eq. Q. Other 
parameters are = 10, CT 2 = 0, 7 = 1 and K = 100. The 
black line indicates the parabola s = j'i, which is our pre¬ 
diction for Pflx = 0.5. The inset shows cuts for exemplary 
values of s = {0, 0.5, 2, 5} in {red, violet, blue, green}. 

Even though environmental variability causes a drift 
term favoring the traits which is less sensitive to envi¬ 
ronmental changes [33], the interplay between drift and 
diffusion term has to be understood to predict the evo¬ 
lutionary outcome. This is even more important as for 
the particular situation discussed here, the environmen¬ 
tal contribution to the drift caused by os is intrinsically 
connected to the diffusion term. Therefore we study the 
fixation probability, i.e. the probability that trait 1 fix¬ 
ates or trait 2 goes extinct. This quantity can be cal¬ 
culated by solving the backward FPE, 0 = C\^Pf,ffxQ) 
for the boundary conditions Pfix(O) = 0 and Paffl) = 1. 
The solutions is given by (for details see SM), 

1 —exp{C[Tanh“^a-|-Tanh“^a(2a; — 1)] } 

1 — exp {2C Tanh“^a} ’ 


with /3 = ^K{ol — 2eoi02 + o^h, a = fij 
and C = 2K (ctJ — o\ — 2s) /(/Jy-^/S + fff). In Fig. 0 we 
show the fixation probability for different values of s and 
CTi = A (ct 2 = 0). Results are obtained by the numerical 
solution of Eq. Q, i.e. before marginalization on A. The 
parabola s = A^2 (Fig. black line) defined by x* = 0.5 
(or Pfix = 0.5), separates the regions where one of the 
two traits is predominant: in the grey (green) area, the 
smaller variability (growth rate) dominates, respectively. 
The general case of both species having variable birth 
rates yields analogous results: a selection advantage for 
the trait with less variability. In the inset, the fixation 
probability depending on A is compared to the analytic 
solution (Eq. (S 6 )) for the four values s = (0,0.5, 2, 5}. 
Both plots demonstrate the advantage of the less variable 
trait. For strong environmental variations it is then ben¬ 
eficial for a species to minimize its sensitivity to those 
variations rather than optimizing its growth rate. In¬ 
terestingly one can interpret this result in the context 
of game theory: Decreasing the sensitivity to environ- 
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mental changes also means to optimize the worst case 
scenario outcome: The average birth rate is the least re¬ 
duced when the variability is small. In game theory, this 
corresponds to the MaxiMin strategy which was shown 
to be very successful in many fields as finance, economy 
or behavioral psychology [SH |3S] . In the context of evo¬ 
lutionary dynamics another example of a MaxiMin strat¬ 
egy was proposed for bacterial chemotaxis where bacte¬ 
ria move move so as to optimize their minimal uptake of 
chemoattractants |36) . 

Besides contributing to the fitness, environmental vari¬ 
ability also influences fixation probability and time in the 
case of neutral evolution, i.e. vi = V 2 and cri = (72 = cr. 
Such analysis is of great interest, as evolution is often 
studied by investigating how neutral mutations evolve 
over time. In recent years fast-sequencing techniques 
made huge amounts of data available, see, e.g., EH, 
which is now analyzed and interpreted by comparison 
to evolutionary models as the Moran or Fisher-Wright 
models [SB]- While the correlation parameter does not 
qualitatively influence results discussed so far, it plays 
an important role for neutral evolution. For fully corre¬ 
lated noise, e = 1, Eq. (S5) is the same as for (ji = (T 2 = 0 
and thereby correspond to the ones obtained for no envi¬ 
ronmental noise, extinctions are solely driven by demo¬ 
graphic fluctuations and well-known results apply [S] . In 
contrast, for all other values of e, including uncoupled 
noise e = 0 , the dynamics differs in two major respects. 
First, the drift term —(1 — e) dx{^ — 2x)x{l — x)P{x) 
does not vanish and corresponds to a stable fixed point 
at X* = 0.5. Second, the diffusion term consists of 
demographic ^x{l — x) and environmental fluctuations 
(1 — e) (T^x^(l — x)^. As the drift suppresses extinction 
events while a larger diffusion term favors them, a more 
detailed analysis is required to grasp the evolutionary 
outcome. 

Due to the stable fixed point, the fixation proba¬ 
bility qualitatively differs from the linear dependence 
= xo which holds for constant or uncorrelated en¬ 
vironments. In Fig[^), typical solutions for e = 1 and 
e = 0 are shown which clearly demonstrated the ensuing 
s-shape for the uncorrelated case (e = 0). Another im¬ 
portant quantity, the extinction time, T(xo), also obeys a 
backward FPE, —1 = £l^T(xo). Employing the bound¬ 
ary conditions r(0) = T(l) = 0, the fixation time in the 
neutral case (s = 0 and tri = (72 = cr) can be calcu¬ 
lated BOl : 
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FIG. 2: Fixation probability (panel a) and time (panel b) 
in the neutral case. Solid lines indicate analytical results for 
the two typical cases of perfectly correlated and uncorrelated 
noises e = 1 and e = 0. Parameters are: = z /2 = 10, 

K — 100, CTi = (72 = 0.5 and 7 = 1 . Dots are simulations of 
the IBM. Additional parameters are m = 1, <j)i = (j )2 = W, 
LOi=uJ 2 = 5, T = 0.01, {E) = 0, ai = Q !2 = 1, and VarfA] = 100. 


nario, = —Ar/ 7 [xo In(xo) + (1 — xq) ln(l — xq)] (see 
Fig#). Fluctuating environments decrease the fixation 
times for all initial conditions. This has a crucial conse¬ 
quence: when measuring extinction times and comparing 
them to standard models without environmental fluctu¬ 
ations, one can only explain large diffusion constants by 
small population sizes. Indeed, it is often found that 
effective population sizes are much smaller than the cen¬ 
sus population sizes [33] ■ Fig. shows that conspic¬ 
uous orders-of-magnitude reductions in the population 
size set in already at moderate levels of environmental 
noise. Amongst other explanations this could account 
for a difference between effective and census population 
size. In other words as long as the level of environmental 
noise and the correlation level of its influence on different 
growth rates is not known, the effective population size 
can only be interpreted as a lower bound for the census 
population size. 

To further investigate the impact of variable environ¬ 
mental conditions, we introduce an exemplary Individual 



T = 


1 


, l-F+(l-Xo) 
.^"i-r_(i-xo) 


ln(l-xo)-l-ln^—In xq 

1 i — Xq 


+Fr+(xo) - Fr_(xo)] , 


(5) 


where a = {l — e)a, C = ^1 -|- F± =2/(l±C'), 

the function Fr(x) = Li 2 (r(l —x)) -b Li 2 (Fx) — Li 2 (r) 
and Li„ is the poly logarithm. The result for e < 1 dif¬ 
fers again from the non-fluctuating/fully correlated sce- 


FIG. 3: Reduction of the effective population size due to en¬ 
vironmental noise in the neutral case for xo = 0.5. Using fsf, 
we plot the values of the population size K and noise a that 
lead to the fixation time T = 100 for 7 = 1 . The black line 
corresponds to (t = 0 (see f5f), i.e. either perfectly correlated 
noise or no environmental noise. In the presence of environ¬ 
mental noise, the values of K are systematically higher and 
increase several orders of magnitude even for moderate noise 
levels. 
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Based Model (IBM). In particular, the IBM serves as a 
proof of principle that linear multiplicative noise can be 
realistically expected and enables us to study the effect 
of such noise beyond the white noise regime. Impor¬ 
tantly, our results presented above hold for any micro¬ 
scopic model whose macroscopic representation is given 
by Eq. Q, i.e. where noise in the birth or death rates 
is linearly multiplicative. In specific scenario discussed 
here, the reproduction rate of an individual, z, at time t, 
depends a priori on the history of environmental condi¬ 
tions experienced during its lifetime where 

fp is the time of birth. This could for example account 
for the level of nutrients or detrimental substances that 
individuals are exposed to. Following [a ED], our model 
is based on independent birth and death rates, now de¬ 
pending on the environmental variations subsumed in the 
scalar value E. The number of environments experienced 
by an individual, z, is denoted as M* and their values are 
contained in a vector E* = (E\,E\, For the 

fully correlated case (e = 1), E is the same for both traits 
while two different values are drawn for e = 0. Environ¬ 
mental conditions change stochastically at rate 1 /r and 
are distributed according to a distribution, p{E), with 
mean {E) and variance Var[E]. 

We first consider a constant environment E. The aver¬ 
age instantaneous growth rate XsiE) is assumed to be a 
positive, monotonically increasing function of E |42| . In 
particular, we consider the sigmoidal function: 

Xs{E) = (j)s+ujs tanh {asE/2) , (6) 

with (ps the ordinate of the inflection point, ujs < (pS the 
maximal deviation from it, and as scales the sensitivity 
to environmental changes. 

Let us now consider changing environments and indi¬ 
viduals whose current growth rate memorizes previously 
experienced environments. The reproduction rate FJ^p,. g 
of an individual, z, of type S now depends on the whole 
vector, E*. For concreteness, we assume that the rate is 

_ M' 

f = l 

where the memory parameter m S [ 0 ; 1 ] defines the in¬ 
fluence of previously experienced environments upon an 
individual’s growth rate. For m = 0 only the current 
environment sets FJ^pj. g = X{E'‘j^i), while for m —)■ 1 all 
experienced environments, M*, have the same influence 
in the arithmetic mean F(gp,,_g = ^ Efc<M> In¬ 

dependent of m we assume that offsprings lose memory 
at the time of reproduction. Bounded growth is modeled 
by death rates g = yiV/iF. 

To compare the results of the microscopic individual 
based model to the effective stochastic model, Eq. Q, 
the parameters of both models have to be mapped. For 
simplicity let us consider the case {E) = 0 and a symmet¬ 
ric distribution p{E) throughout the following discussion. 



0.01 1 100 10000 

environmental rate, 1 /t 

FIG. 4: Comparison of the IBM and the Langevin model. We 
show the fixation time for neutral evolution for xo = 0.5 vs 
the environmental switching rate 1/r. Dots correspond to the 
IBM [m = 0] for different values of cui = tU 2 = {0.5, 2, 5} in 
red, blue and green. Black lines are analytic solutions [Eq. ([^ 
with Eq. Qj. For quickly fluctuating environments both re¬ 
sults are in good agreement whilst for large r the white noise 
approximation fails. Other parameters are as in Fig. 


Since death rates are constant, there is a direct corre¬ 
spondence between them in the Langevin and the IBM. 
For birth rates and their STDs the situation is more in¬ 
tricate as we discuss hereafter. For the no-memory case 
(to = 0) an exact mapping is obtained [SD]: For strong 
fluctuations, a|Var[E] ^ 1, the mean of the growth rate 
vs and STD of the noise as in Eq. Q are given by: 

vs{m = 0) = (/)5 -I- cts(to = 0) = wV^. (8) 


Note that the variability in the growth rate not only re¬ 
sults in as > 0 , but also influences the average repro¬ 
duction rate us- While for to = 0 such a variability 
increases vsi the second term of vs is reduced while to 
increases till it changes sign (see SM [30] for details). 
For instance, for to = 1 the growth rate is approximately 
(pS ~ wP'T. Hence, the more variable trait has a disadvan¬ 
tage in the average reproduction rate in addition to the 
effects discussed above. For to = 1, the approximation 
as{m = 1) « w!s\/t holds |30|. Dependencies in this ex¬ 
pression are intuited as follows. The number of environ¬ 
mental changes an individual experiences until the mem¬ 
ory resets is of the order M ~ tufe/r, where tufe oc l/z /5 is 
the typical time for an individual to reproduce or die. As 
environmental changes are independent random events, 
the variance of the reproduction rates (S21) is oc ujg/M. 
The expression for as{m = 1) is finally obtained noting 
that correlations in the noise extend over times ~ tufe 
therefore it follows that the average reproduction rate vs 
drops out. 

For a detailed comparison of the IBM with the analyt¬ 
ics derived in the first part of this paper, we simulate the 
IBM with a modified Gillespie algorithm updating repro¬ 
duction rates after every environmental change |27j . As 
shown in Figs.j^) and b), results for fixation probability 
and time, are in very good agreement with analytic so¬ 
lutions [Eqs. ( |S^ and ([^]. In particular, the sigmoidal 
shape of the fixation probability is well reproduced by the 
IBM, supporting the existence and importance of linear 
multiplicative noise. 









5 


Finally, the IBM enables us to study the environmen¬ 
tal switching rate. This is of main interest as results 
obtained previously strictly only hold for very rapidly 
fluctuating environments. In Fig. the dependency on 
r of the extinction time in the neutral case for xq = 0.5 is 
shown for different ojs; see SM for results with s ^ 0 [50] . 
The black lines correspond to Eq. ([^ mapped according 
to Eqs. ^ and dots are obtained by stochastic simula¬ 
tions of the IBM for m = 0. For t < 1 both models 
are in very good agreement. This demonstrates that the 
white noise approximation is valid in a broad parame¬ 
ter range, where fluctuating environments substantially 
influence the evolutionary dynamics. 

In summary, we demonstrated that environmental 
variability has crucial impact on evolutionary fitness. 
First, we quantified the role of reduced sensitivity to 
environmental changes and determined how it substan¬ 
tially increases the fitness. Second, we showed that the 
timescale of extinction in the neutral case is strongly af¬ 
fected by environmental noise. That provides a mecha¬ 
nism to explain experimental observations of population 
sizes that are often much smaller than expected. Finally, 
we investigated individual based models that generate 
the linear multiplicative noise considered here. It will be 
of interest to investigate how different forms of memory 
or time-dependent reproduction rates influence evolution 
and to integrate them with evolutionary game models. 

We thank Jonas Cremer for valuable discussions and 
comments on the manuscript. AM acknowledges the Ger¬ 
man Academic Exchange Service (DAAD) for financial 
support. 
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Supporting Information 

In this Supporting Material, we provide more de¬ 
tails on the calculations leading to the one-dimensional 
Fokker-Planck equation (FPE), Eq. (3) main text. 
Moreover, we present calculations for the average 
fixation time and the extinction probability. We discuss 
the mapping of the individual based model (IBM) to 
the Langevin model. For the no-memory limit we 
present an analytic derivation of the mapping. For other 
parameter values we give heuristic arguments that are 
supported by additional data. Finally, we present results 
for non-neutral evolution investigating the regime in 
which the white noise approximation is an adequate 
description for the evolutionary process. 


DERIVATION OF THE ONE-DIMENSIONAL 
FOKKER-PLANCK EQUATION 


In this Section, we provide details how the Langevin 
equations, 

Ni = - 7^) + y^^vr(iq"TT^/.ti-l-criIVi^i 

N 2 = N2{v2 - 1%) + \JN2{v2 + 7^ )M 2 +CT 2 iV 26 (Sl) 

can be transformed into the one-dimensional FPE pre¬ 
sented in the main text. From general results on 
stochastic processes (see [ 1 ]), it follows that the previ¬ 
ous Langevin equation is associated to the following two- 
dimensional FPE: 


5P(7Vi,V2,t) 

dt 


= edl 2(^ia2NiN2P - ^ 




p 


N 


(S2) 


where di = 8^.. The drift part is directly stemming 
from the non-fluctuating parts of the Langevin equations 
Ns{vs — "fN/K). Diffusion depends on the correlation 
level of the noises experienced by the two species. In 
particular, we have introduced the correlation coefficient 
e = (Ci? 2 )/\/(Ci)(^ 2 )' The case when the two noises are 
the same is given by e = I, when they are independent is 
e = 0 and when they are anti-correlated is e = — I. 

To study the evolutionary dynamics associated to 


Eq. (S2), the relative abundances are the natural choice 


of variables. Therefore, we transform the absolute abun¬ 
dances Ni and N 2 to x = and N = Ni + N 2 . 

To perform the change of variables, not only Ni = xN 
and N 2 = (1 — x)N have to be replaced, also the dif¬ 
ferential operators and the probability distribution have 
to be transformed. Ensuring that the latter is still nor¬ 
malized after change of variables, the Jacobian has to be 
introduced, P(iVi, V 2 ) —)■ j^P{x, N). The derivatives are 
given by, Oni ^ + On and dN 2 ^ -§dj; + On- 


After the change of variables, the FPE for x and N can 
now be further simplified exploiting the fact that the time 
scale of selection, s = 121 — 122 , is much slower than the one 
of the population growth z 2 ix+:^ 2 ( 1 — 2 ;) [2]. Therefore, we 
marginalize the FPE with respect to the total population 
size N. Thereby, the integrals dN of V-derivative 
terms such as On* or Nd%* = dn {NOn*) — On* vanish 
and the FPE simplifies to 


dP{x, t) 


= dx{ [—s — cr|( 1 —a;)-I--I-ecTi (72 (1 — 2 a:)] Q} 


+a 


dt 

&) 

ctJ — 2 ecritT 2 +o'| 


n ^ , 7 , - sa: 


Q 


(S3) 


where Q = x(l — x)P{x,t). The drift term in the sec¬ 
ond line stemming from demographic fluctuations can be 
neglected as ^ I holds. To finally arrive at the one¬ 
dimensional FPE employed in the main text, we compute 
the steady state population size N*. As the deterministic 
differential equation for N is given by 


N = N 


X 12 l -I- (1 - X)I 22 - 7 — 


the fixed point for the populations size is N* = 
K/j[i 2 ix + 122(1 — x)]. Employing that relation and the 
aforementioned condition s <C vix + 122(1 — x), the last 
term in Eq. ( [S^ can be simplified as ~ 2 K> which 

finally leads to the one-dimensional FPE in the main text: 


dP(x, t) 


dt 


= dx{ \^—s — a2(l — x)+a\x + e(Jia2 (1 — 2x)] Q} 


+dt 




Q 


(S4) 


FIXATION PROBABILITY 


In the following, we derive a general expression for the 
fixation probability. The calculations are analogous to 
the procedure for the neutral case described in the body 
of the paper. To determine the fixation probability the 
following backward equation has to be solved, 


0 =x (I 



- x) I [s 


|S + (T2(l—x) — Cr^X — 
— 2 eCTi (72-fcr? 


-x(I-x) 


e(7i(72(l-2x)] dx 

a2|pfi,,(x). (S5) 


Boundary conditions are Pfix(O) = 0 and Pfix(l) = 1- 
The solution to Eq. (|S5[) for the fixation probability is 


J^fix(x) 


I —exp{C[Tanh ^a-l-Tanh ^a(2x—1)]} 

1 — exp { 2 (^Tanh“^Q!} ’ 

(S 6 ) 
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with 

/? 


change of variables) or first integrating by parts to sat- 


c = 


jK{al-2t(Jia2 + (rl)h\ a=—^===- 

y/8 + /3^ 

2K (cr^ — — 2s) 

/37\/8TF^‘ 


isfy the condition a > 0 in (Sill. The resulting ex- 


(S7) 


The solution (S 6 ) is obtained by integrating (S5) once, 


to find the gradient 

dMx) = const. + . (S 8 ) 

(1 - a(2a; - 1))^/"+^ 


pression is then transformed to the form given in the 
main text (which is the one given by Mathematica) 
by using the reflection property, Li 2 (a:) + Li 2 (l — x) = 
Li 2 (l) — lnxln(l — x), see [5]. 


General case 

In the general case when selection is present, the ex¬ 
pression for the average fixation time cannot be found 
explicitly but is reducible to quadratures as follows. The 


The expression (S 81 is verified to be proportional to the fixation time obeys the backward equation (S5) with the 


derivative of ^^ boundary conditions 

are then imposed to fix the two constants of integra¬ 
tion. The resulting expression is finally transformed into 
Eq. ( [S^ by using the elementary identity: 2 Tanh“^ (x) = 
log [(1 -I- a;)/(l — a:)]. It is verified that in the limit C —t 0, 
one recovers the expression given in the main text. 

All in all, the behavior we discussed in the main text 
is validated by analyzing the fixation probability: Both a 
higher growth rate and a smaller variability are beneficial 
for an individual. 


AVERAGE TIME FOR FIXATION 
Neutral case 

The expression for the time of fixation in the neutral 
case that we presented in the body of the paper is de¬ 
rived as follows. The average time for fixation obeys the 
following backward equation, 

{l-2a:-k [x{l-x) d^} da:T{x) = 

— (a^x{l — x)) (S9) 

with d = (1 —e)cr the boundary conditions T(0) = T(l) = 


0. Integrating Eq. (S91 and by variation of constants, we 
obtain: 


dxT{x) = 


1 




1 f 1 — X 
A -I- In 

/T^ 


,(S10) 


where A is a constant to be fixed by the boundary condi¬ 


tions. The integrals of Eq. (SIO) needed for T(x) are 
performed by decomposing the rational function at the 
prefactor and using the formula: 


In (a -I- 6a;) , , , t ■ ^ 

-da: = In a In a: — L 12 - 

X \ a 


a > 0 , 


(sii) 

that follows from the very definition of the dilogarithm 
Li 2 (a:) = — ln{l — u)/udu (see [3]). The formula 
(Sll) is used four times either directly (with a simple 


left-hand side replaced by —I. Using the definitions (S71, 
we obtain 

■|2 —(s -I- 0-2 — e(TiCT 2 ) — 2/3^x -I- [/3^a:(l —x) -I- 2] x 


dxT{x) = - 


2K 


7 x(l — x) 


(S12) 


Boundary conditions are T{Q) = r(l) = 0. The homo¬ 
geneous solution was already found following (S 8 I and 
reads 


Thom{x) = Cl +C2 


( X+{x) 


C/2 


(S13) 


\X-{x) 

where Ci and C 2 are constants and we defined 

X+{x) = l + a{2x-l), x_(x) = l-a(2x-1), (SI4) 

to simplify notation. The non-homogeneous solution for 
the gradient of T is obtained by varying the constant in 


(S 8 ), remarking that j3'^x{l — x) + 2 = x+(x)x_(x )(8 -I- 
and integrating the resulting first-order differential 
equation to obtain 


^xd^partix^ — 


A:x-^(x)^/^ ^ 
2C/27(C/2 + 1) 


-k 


(a - l)Fi + 1 , 1 ;^ + 2 ; ^ 

(a + l)Fi (^ + ^ 


where Fi is the hypergeometric function of two vari¬ 
ables [1]. The solution for T involves the integral 
lo^dy Tpart (y) dy of the expression above (for which a 
closed form does not seem to be available), and the two 
constants in (SIS) are fixed by 


Cl + C 2 dy. 



































on E is given by the sigmoidal function: 



std, A 


FIG. SI: Extinction time depending on A for different values 
of the selection strength: s = 0.5 (red), s = 2 (violet) and 
s = 5 (blue). Dots are numerical solutions of the Langevin 
equations, Eq. (SI I, and black lines are solutions of Eq. (S12 1 . 


As(£^) = + wstanh , (S15) 

which reduces to (f>s i ujs in the limit of large variances 
Var[E]. Birth rates are dehned as, 

_ M' 

r:.p,,s = (S16) 

In the no-memory limit m = 0, the growth rate is there¬ 
fore given by the instantaneous growth rate Xs{E), while 
for TO —>■ 1 the current growth rate is the arithmetic mean 
of all previously experienced environments. Death rates 
are given by T^ae!,th,s = l^/K. 


It is verified from the expression above or directly from 
the original equation ( S12[ ) that in the two limits x —>■ 0 
and a: —?> 1 the solution behaves like in the neutral case, 
i.e. —K/jx\ogx and — a;) log(l — a:). Selec¬ 

tion and the rest of the parameters affect of course the 
solution in the rest of the interval of definition x G [0,1]. 


COEXISTENCE TIME 


Depending on the position of the stable fixed point, co¬ 
existence between two species (one with a larger growth 
rate, one with a smaller variability, i^i > V 2 and tJi = 
A, (72 = 0) is possible. In this section we present some 
additional data demonstrating this. In Fig. m the ex¬ 
tinction time which corresponds to the time of coexis¬ 
tence is shown depending on A is shown for different 
values of s. Dots correspond to solutions of Eqs. (SI) 
and black lines are numerical solutions of Eq. (S12). The 
extinction time has a maximum which exactly coincides 
with the parameter values of a hxed point x* = 0.5. The 
dependence of this maximal extinction time on the selec¬ 


tion strength s is shown in Fig. S2 


MAPPING INDIVIDUAL-BASED MODELS 
ONTO THE LANGEVIN DYNAMICS 


The aim of this Section is to show that individual- 
based models are described by the Langevin equations. 


Eqs. (SI), discussed in the main text and to analyze the 


mapping between the parameters of the two models. 

The environmental conditions change stochastically at 
the rate 1/r and are distributed according to a distribu¬ 
tion, p{E), with mean (E) and variance Var[E]. The de¬ 
pendency of the instantaneous reproduction rate Xs{E) 


No memory, m = 0 


We discuss first the model without memory, where the 
memory parameter, to, is zero: Individuals reproduce 


with the instantaneous reproduction rates [Eq. (S15)], 
which reduce to (psEtos in the limit of large environmen¬ 
tal variance. We consider an interval of length 5t t 
such that the probability for an individual to reproduce 
or die is small, yet the total number of events occurring 
over the whole population (^ A 3> 1) is large. Neglecting 
the standard demographic noise term [5] , the variation of 
the ^-type population is given by 

Ns{t + 5t) ~ Ns{t) + Nsit) - 7 ®) 5t+ 

pt-\-St 

-l-ws / Ns{s)a{s) ds 


(S17) 



FIG. S2: Extinction time for different values of s and A = 
■\/(2s). This combination of s and A corresponds to a stable 
fixed point at x* = 0.5 and the maximal coexistence time for 
each value of s, see Fig. |Sl| As not only the selection strength 
but also the variability is increasing from left to right, the 
fixation time is a monotonically decreasing function of s. 
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where (t(s) is the environmental Boolean random vari¬ 
able that takes values ±1 and switches with character¬ 


istic time T. The last term of Eq. (S17) is estimated as 
follows 

pSt 

/ Ns{t + s)a{t + s) (is ~ Ns{t)G^-\- 

Jo 


p6t ps 

+Ns{s)uJs / a{t + s)ds / d{t + s')ds', 
Jo Jo 


(S18) 


where is a Gaussian random variable having zero mean 
and variance 


n6t p6t 

Var[0^] = / ds ds'{d-{s)d{s')) = . 

Jo Jo 


(S19) 


Here, we used that {a{t)a{t')) = e I* * 1/"^ and St ^ 
T. The second term in Eq. (S18) is evaluated at the 


order St using the same integral, Eq. (S19), and gives 


Ns{t)u}sT6t. Combining back all the terms, we conclude 


that the equation (S17) is equivalent to the Langevin 
equation ( |S1[ ) with the mapping of the parameters 

>^3 = ( I)S + w|t ; ag = 2w|t . (S20) 

Note that the standard demographic noise term in 


Eq. (SI) should a priori include the fluctuating environ¬ 


mental term Ns{t)ujsQ‘^ in the sum of the rates. In fact, 
it can be safely ignored as (JigNsSt ^ ijjgNg^/2T5t due 
to 4>g > ujg and 5t ^ r. 


Einally, the factor 2 appearing in ct| in (S20) depends 
on the Poisson statistics of the environmental fluctua¬ 


tions. If the duration is fixed and equal to r, Eq. (S19) 


becomes rSt. In that case, the corresponding mappings 
are vs = 4>s -\- LOgT 12 and (t| = w|r. This is confirmed 
numerically in Fig. |S3| where we show data for exponen¬ 
tially distributed (black) and fixed duration (red) envi¬ 
ronments. Solid lines are analytic solutions of the fixa¬ 
tion time [Eq. (6) main text] employing the respective 
mappings. 


Finite memory, m > 0 

We now turn to the scenario where memory extends 
over several environmental conditions that an individual 
previously experienced. Whilst for m = 0 an exact an¬ 
alytic mapping can be found in the limit of small r, for 
finite memories the situation is more intricate. However, 
for the special case to = 1 the variability in the growth 
rate can be well approximated by the following argument. 
The reproduction rate TJ^pj. g of an individual, i, of type 
S at time t depends on the average of all instantaneous 
reproduction rates Xg{E) experienced by the individual: 

rLpr.s(i) = i E ^siEl) = ^s+ TW) ■ (S21) 



FIG. S3: Comparison of data obtained by simulations in the 
neutral case with fixed (red) and exponentially distributed 
(black) environmental changes. Both sets of data agree with 
our analytic calculations, where we used the mappings = 
Lo'^T for fixed times of environmental changes and = lo‘^2t 
for exponentially distributed switches. Thereby, the data con¬ 
firms that the origin of the factor 2 in the mapping is solely 
the exponential distribution of the environmental changes. 
Parameters are 01 = 02 = 1 and Wi = 0^2 = 0.9, 7 = 1 , 
K = 5000, {E) = 0, VarlF] = 100 and r = 0.01. 


At the time of reproduction, we assume for simplicity 
that offspring looses its memory of past environments 
experienced by the progenitor. 

We consider again a time interval of length St as in 


(S17). Fluctuations in the rates Eg{t) decorrelate on 
timescales of the order of the lifetimes of individuals, 
tiife, which are much longer than r and St. Therefore 


on the St scale, noise is smooth, contrary to (S17|. Con¬ 


versely, timescales of several lifetimes are much smaller 
than those on which selection acts and much longer than 
the characteristic time of the noise. Therefore, to de¬ 
scribe the dynamics of the fractions, the environmen¬ 
tal noises are well approximated by a shortly correlated 
noise. An estimation of the amplitude of the noise is 
obtained by calculating the sum 




1 


00 Ns Ns 


- E EE(rs(4)A4f^g(t,)At,). 


.^= — 00 2=1 j — 1 


The durations At of the environmental intervals are inde¬ 
pendent for different k and £ (and the contribution k = £ 
is negligible with respect to the rest of the sum) so that 
one can replace them by r. In addition, the symmetry in 
the indices of the intervals allows us to further simplify 
the expression 


<^S 


00 Ns Ns 

2^EEE(rM^fc)f'5(^^))- (S22) 


£—k i—1 j—1 


To compute the average in Eq. (S22) three different or¬ 


ders of events have to be distinguished: a) If the birth of 
the j-th individual was prior to the one of the i-th indi¬ 


vidual, then it follows from Eq (S21) that the quantity to 
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be averaged is w|M^/ where Mj is 

the number of environmental switches since the birth of 


the t-th individual up to time tfc. To derive Eq. (S221 we 


have used that the terms in the sum (S21) take indepen¬ 


dent values iws with equal probability. Conversely, case 
(b) is when the birth of the j individual is posterior to the 
one of the Tth individual. Then the quantity to be av¬ 
eraged is LUg {^Ml — 6b) / where 5b is the time 

between the birth of the z-th and the j-th individuals. 
Finally, in case (c) when 5b > Ml, the correlation is zero 
as there is no overlap between the environmental fluctua¬ 
tions of the two individuals. Due to the Poissonian nature 
of the events, the number of switches since birth (back 
in the past) or before reproduction (forward in the fu¬ 
ture) have the same distribution exp{—M/M)/M where 
M ~ tiife/r (its exact value does not affect the sequel). 
It follows that 


2 2a;|r 

M 


POO poo 


1 


, dt du 
lo Jo Jo u + v + t 

^ u — V 


+ dt du 

Jo Jo Jo [u-v + tju 

(S23) 

The integral over t is the continuous approximation of 


the sum over i — k appearing in Eq. (S22) while the vari¬ 


ables u and V refer to the variables Ml and 5b. The first 


and second term in the square parentheses of (S231 cor¬ 


respond to cases (a) and (b), respectively. By a series of 
change of variables and integrations by parts, it is shown 
that (S231 reduces to 




2u}i 


M 


du 


^-u/M 


dv e = iVgT . 


(S24) 


The validity of this approximation is confirmed for the 
neutral case in Fig. 2 in the main body of the paper where 
the fixation probability and time are compared to the 
analytic calculations employing Eq. (S24). Additional 


data for non-neutral evolution is presented in Fig. S4 
where two species (one with finite variability, uj = 0.9, 
one with vanishing variability) are analyzed. Analytic 
solutions for the fixation time and probability are fitted 
to simulation data. The best fit deviates less than 1% 


from Eq. (S24|. 


We now briefly discuss the dependence of as on the 


memory parameter m. In Fig. S5 the STD of the noise in 
the growth rate, a ,depending on the memory parameter, 
m, is analyzed. Results were obtained by simulating the 
neutral evolution case, measuring the fixation time and 
calculating cr employing the analytic expression for the 
extinction time [Eq. 6 main text]. For to = 0 the thereby 
obtained value agrees nicely with the mapping introduced 


above indicated by the red dashed line [Eq. (S20)]. With 
increasing memory, to, the STD of the noise, a, decreases 
monotonically and approaches Eq. ( S24[) for to —> 1. 



FIG. S4: Fixation probability and time for non-neutral evo¬ 
lution with memory m = 1. Black dots correspond to simu¬ 
lation results of the IBM and red lines are analytic solutions. 
To obtain the latter we fitted Eq. (S 61 and the solution of 
Eq. (S12| to the IBM. We used A — s + — cria 2 e and 

B = af + 2 (Ti(T 2 e + nl as fitting parameters and obtained 


A = -3.04 X 10-^ and B = 8.25 x 10” 


Other parameters 


are (l>i = (l >2 ~ 1, uti = 0.9, a ;2 = 0, r = 1/100, 7 = 1, 
K = 5000, (E) = 0 and Var[i5]. 


Let us now analyze the mapping of the average re¬ 
production rate vs- Importantly, this mapping is very 
sensitive to model details which we will exemplify in the 
following. As results for neutral species do not depen¬ 
dent on the average reproduction rate, we have to turn 
to the evolution of non-equal individuals to understand 
the mapping of the growth rates. In Fig. we show 
the fixation probability for two species with the same 
4>i = 4>2 = ^ but only the first species has a variable 
reproduction rate uji = 0.9, ClI 2 = 0 for to = 1. Red 
lines correspond to a fit with s = z/i — 1/2 ~ —0.0030 and 
agree perfectly with simulation results. In other words, 
the first species does not only have a disadvantage due its 
sensitivity on environmental changes, ai > a 2 , but also 
has a smaller average growth rate. To study this effect 


b 


0.3 


(1 __ 


0.25 \-° o . 
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memory, m 


FIG. S5: Numerical estimation of the STD of the noise depen¬ 
dence on the memory parameter m. Extinction times in the 
neutral scenario were measured. By evaluating the inverse of 
the fixation time function [Eq. ( 6 ) main text] the STD of the 
growth was calculated. For m = 0 the result agrees with the 
mapping, cf. Eq (S20l, indicated by the red dashed line. For 


larger values of m the variability is reduced. Parameters are 
</>i = 02 = 10, tJi = tJ 2 = 8.237, r = 1/500, 7 = 1, A = 100, 
(E) = 0 and Var[A]. 
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0.2 0.4 0.6 

memory, m 


(b) 



FIG. S 6 : (a) Fixation probability and (b) selection coefficient 
depending on the memory parameter m. The first species’ 
growth rate depends on the environment while the second 
one’s is constant. Black dots correspond to the IBM intro¬ 
duced in the main body of the paper, while red dots represent 
a model modification not memorizing the first environment 
(the one in which an individual is born; for details see text). 
While the fixation probability is a direct simulation result, the 
selection coefficient s is inferred from it using the additional 
data presented in Fig. |S5| Parameters are 4>i = 4^2 = 10, 
a;i = 8.237, tuz = 0, t = 1/500, 7 = 1 , if = 100, {E) = 0 and 
Var[i;]. 


in more detail, let us now analyze the fixation probabil¬ 
ity dependence on the memory parameter m, see Fig. |S 6 | 
panel (a). Black dots correspond to the standard IBM 
(if not mentioned otherwise our discussion applies to this 
data), red dots to a slightly changed model which is go¬ 
ing to be introduced in the following. For m = 0 both 
species are equally likely to fixate as the growth advan¬ 
tage of the more variable species i^i = c/) + uj^t > V 2 = (j) 
exactly compensates for its disadvantage due to the STD 
of the noise cti > <72- For increasing values of m first the 
more variable (m < 0.7) later the less variable species is 
favored (m > 0.7). Whether this behavior is caused by 
the STD of the noise or differences in the mean reproduc¬ 
tion rates is not obvious as the influence of both fitness 
contributions is of comparable strength. Therefore, we 
estimate the selection coefficient, s = vi — 122 , from the 
fixation probability data, see Fig. S 6 panel (b). This is 
achieved by assuming that the variability of species 1 
with oji = 0 is zero (tri = 0 ) and that the variability of 
species 2 with uj 2 = 8.237 is the same as in the neutral 
evolution scenario and thereby given by the data pre¬ 
sented in Fig |S5[ Note that this approximation might 
neglect some higher noise correlations arising due to the 
coupling of both species via the carrying capacity. For 
m = 0 the thereby obtained value of s agrees well with 


our analytic results, cf. Eq. (S20). With increasing m 
the growth rate of the more variable species is decreased 
till the selection coefficients becomes negative effectively 
favoring the more variable species. However the decrease 
of the selection coefficient with m is smaller than the re¬ 
duction of a shown in Fig. S5 Hence, for small m the 
more variable species is favored as its advantage due to a 
larger average reproduction rate is larger than its disad¬ 
vantage due to its sensitivity on the environment. This 
advantage in the growth rate is more sensitive to details 
in the IBM in comparison to the variability discussed in 
the main text for the Langevin equation Eq. (l).We are 
going to illustrate this by analyzing a slightly modified 
version of the IBM. But before doing so, we present an 
intuitive argument explaining one factor influencing the 
average growth rate: When an individual is born it expe¬ 
riences the current environment shorter than the average 
length of an environment. However, the model weights all 
experienced environments equally, see Eq. (7) in the main 
text. As the first experienced environment is more likely 
to be a good environment [more reproduction events hap¬ 
pen during more beneficial environments], higher growth 
rates have a larger weight in the average and the aver¬ 
age growth rate of the variable species is effectively in¬ 
creased. To obtain a description including this factor, it 
would be best to perform a time average over all previ¬ 
ously experienced environments. Unfortunately, such a 
procedure is computationally very expensive. We there¬ 
fore, test our explanation for the bias by not including 
the very hrst environment, the one in which an individual 
is born, in the averaged reproduction rate. The red dots 
in Fig. |S 6 | corresponds to simulation results for this mod¬ 
ified model. Even though all parameters are the same 
and for r = 1/500 and c/i « 10 an individual experiences 
in average 25 different environments, the small modifi¬ 
cation of the model substantially changes the simulation 
results. While the modification almost has no impact 
on the STD of the noise, it alters the average reproduc¬ 
tion rate. For instance the regime in which the more 
variable species is favored completely disappears, cf. red 
dots Fig. |S 6 | where Pfix < 0.5 for all m. This example 
illustrates that on the first sight tiny details of an IBM 
might substantially influence the evolutionary outcome 
and that one should be cautious when drawing conclu¬ 
sions from them. Importantly, the mechanism discussed 
in the main text does not rely on specific assumptions of 
the microscopic models: A finite STD of the growth rate 
is always a disadvantage. It might be compensated for by 
a larger average reproduction rate but the same value of 
the growth rate without variability is always preferable. 


DEPENDENCE ON THE SWITCHING RATE 


In this Section, we present additional data for the non¬ 


neutral case. Fig. S7 shows the fixation probability de- 
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FIG. S7: Dependence of the fixation probability on the 

environmental switching rate. Both species have the same 
01 = <?!>2 = 10 , but only the second species is sensitive to the 
environment (ui =0, 012 = 9). For no memory (m = 0), 
both species are equally likely to fixate, as the advantage in 
the average growth rate of species 2 exactly compensates for 
its disadvantage due to its sensitivity on the environment. 
For m > 0 those two effects do not cancel out anymore and 
the second species is favored. Other parameters are 7 = 1 , 
K = 100, a — 1, {E) — 0 and Var[i?] = 100. 


pending on the environmental switching rate 1/r. In 
particular, we investigate extinction for a species which 
is not sensitive to its environment (0i = 0 = 10, uji = 0) 
competing with a sensitive species (02 = 0 = 10, 012 = 9) 
for different values of m. In the case of no memory m = 0 
both species are equally likely to fixate as the advantage 


in the average reproduction rate 1^2 = 0 + uj 2 T exactly 
compensates for the disadvantage due to the STD of the 
noise in the growth rate [see Eq. (S20)]. For larger values 
of the memory parameter, a bias favoring the species with 
w = 0 is present (the exact value of the fixation proba¬ 
bility depends on mapping details as discussed above). 
Importantly, the bias is not only present for very quickly 
fluctuating environments, but already emerges if repro¬ 
duction events happen on a time scale comparable to r. 
This supports the conclusion, that we were already draw¬ 
ing in the body of this paper when discussing Fig. 4: the 
white noise approximation is an adequate description for 
such an evolutionary process in that parameter regime. 
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